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Program Summary 

Title of program: EXCALIBUR 

Program obtainable from: R. Kleiss, NIKHEF-H, 
RO.B. 41882, 1009 DB Amsterdam, The Netherlands, 
t30@nikhefh.nikhef.nl 

R. Pittau, Instituut-Lorentz, University of Leiden, 
RO.B. 9506, 2300 RA Leiden, The Netherlands, 
rulgmO@leidenuniv.nl 

Programming language used: FORTRAN 77 



Memory required: about 170kbytes 

number of bits per word: 32 

Subprograms used: none 

Number of lines in distributed program: 3784 



Keywords: Decaying vector-boson production, all four fermion processes, 
electroweak and QCD background, initial state QED radiation, multichan- 
nel Monte Carlo approach. 



Nature of physical problem: Heavy vector boson production will be investi- 
gated at e + e" colliders in a wide range of energies. At LEP II, the relevant 
process is 

e+e" -> W + W~ . (1) 
At higher energies other processes like 

e + e- -> ZZ , (2) 

e+e" -> Weu e , (3) 

e+e" -> Ze + e~ , (4) 

e + e" -> Z v e v e , (5) 
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become important. The detected experimental signal for all above processes 
is a four fermion final state. Therefore, a Monte Carlo program being able 
to take into account both signal and background electroweak diagrams for 
all four fermion processes is required. QED initial state radiation and QCD 
background play also an important role and have to be included. 

Method of solution: An event generator is the most suitable choice for a 
program to be able to deal with the above physical problem, since each 
generated event is a complete description of the momenta of the produced 
particles and any experimental cut can be easily implemented. There are 
two basic difficulties. First of all the number of Feynman diagrams can be 
very large. Secondly, taking into account also the background diagrams, 
the peaking structure of the matrix element squared is very rich, so that a 
straightforward integration over the allowed phase space is impractical. The 
former problem can be solved by using spinorial techniques to compute the 
amplitudes and taking massless fermions. The latter requires the use of a 
multichannel approach, where the integration variables are generated accord- 
ing to distributions that approximately reproduce the peaking behaviour of 
the integrand, so reducing the estimated Monte Carlo error. 

Since one wishes to take into account all possible final state (that means 
to have from 3 to 144 different Feynman diagrams, many of them leading to 
different peaks in the phase space), a systematic and automatic procedure for 
both the generation of the Feynman diagrams and the phase space integration 
is unavoidable, together with an algorithm for the self-optimization of the 
predetermined probabilities used to choose the various channels. 

All that has been implemented in EXCALIBUR. This paper serves also as an 
example of the entire procedure to be used to build future event generators. 

Typical running time: about 100 events per second on HP, depending on the 
chosen physical process. 

Unusual features of the program: none 
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Long Write-Up 



1 Introduction 

In the near future LEP II will become operative in the energy region around 
200 GeV. The physics relevant at higher energies will be investigated at the 
next generation of e + e~ linear colliders. Many interesting physics issues can 
be studied and one of them is the gauge-boson production. Around 200 GeV 
events with the signature of two produced Ws have a large cross section, 
while single boson production processes become important with increasing 
energy M. One can distinguish five sizeable reactions (eqs. (|I])-(E|)) in which 
gauge bosons are produced. Due to the fact that the massive bosons are 
unstable particles, all those processes end up with a detectable 4-fermion final 
state to which many Feynman diagrams can contribute. Some of them are 
related to the reactions dU)-© (signal diagrams); others are not (background 
diagrams). For this reason a precise knowledge of all possible processes 

e + e~ — > 4 fermions (6) 

is unavoidable in order to make comparison with experiment 0. 

In addition to these background effects, one wants to be able to study 
any experimental distribution, taking into account the dominant radiative 
corrections effects, and the possibility to implement any experimental cut. 
To solve these problems we wrote an event generator, that can handle all 
diagrams leading to a specified 4-fermion final state (with, of course, the 
option of a restriction to the signal diagrams), and that incorporates the LL 
O (a) and O (a 2 ) initial state radiation (ISR), with exponentiation of the 
remaining soft-photon effects ||. Furthermore, with a four quark final state, 
QCD diagrams are present as well, giving non-negligible effects that have 
been also included [§]]. 

It should be noted that even when cross sections do not dramatically 
change under inclusion of tiny effects, there are quantities that are very sen- 
sitive to any small correction. Among them is the average energy e radiated 
by the beams. A precise knowledge of e is required at LEP II when the 
reconstruction of the jet invariant mass distributions is performed to mea- 
sure the W mass |§, §. In addition, e is also very sensitive to the imposed 
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experimental cuts ||, so that, once more, a Monte Carlo approach is to be 
preferred. 

In order to build a fast program we have taken the limit of vanishing 
fermion masses. Even if this implies the absence of diagrams where a Higgs 
boson couples to the fermions - and therefore we cannot compute the Higgs 
signal - we can at least estimate the background. On the other hand, the 
inclusion of the leading Higgs signal is trivial because only few diagrams 
account for it and, due to their helicity structure, do not interfere with all 
the others in the limit of massless fermions. However, this has not been 
implemented. 

In the rest of this paper we shall describe EXCALIBUR, our event generator 
to compute all 4-fermion processes in e + e~ collisions, including QED initial 
state corrections and QCD diagrams. The general structure of the code is 
flexible enough to deal with physics at the energy scales from 100 GeV to 1 
TeV. 

2 Theory and general features 

There are two sources of complications. First af all one has to generate and 
compute all possible Feynman diagrams contributing to any given final state. 
Then the Monte Carlo integration has to be performed. 

As explained in ref. || the former problem can be efficiently solved by us- 
ing spinorial helicity techniques. The amplitudes receive contributions from 
Abelian and non-Abelian graphs, with two distinct topological structures 
(see fig. 1). In these so-called generic diagrams, all particles are assumed to 
be outgoing: assigning two fermion legs to be the initial-state fermions (by 
crossing), the actual Feynman diagrams are generated. The particles and 
antiparticles can each be assigned in six ways to the external lines (in prin- 
ciple). This gives 36 possible permutations. The Abelian diagrams are built 
by selecting, for each permutation, only those cases in which the exchanged 
bosons, that may be W + , W~, Z or 7, give rise to existing and charge con- 
serving vertices. In the non-Abelian diagrams, two of the vector bosons are 
fixed to be W + and W~, and the third one can be Z or 7. 
This procedure gives, for the Abelian graphs, a maximum of 144 different 
diagrams, and at most 8 for the non-Abelian diagrams. 
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Figure 1: generic diagrams for four-fermion production. The fermion momenta 
and helicities, and the bosons are indicated. The bosons can be either Z, 
W^, or 7; V can be either Z or 7. 

The spinorial structure of each diagram can always be written in such a 
way that a particular combination of axial and vector couplings factorizes 
for a given helicity assignment. For example, if a, and Vi are the axial and 
vector couplings in the vertices of the abelian diagram of fig. 1, the following 
equation holds 

u x {l)Y{vi + a ll5 )u x {2) x (7) 

u P (3)j^(v 2 + 0275) (A +^2+^3)7^3 + a 3 7 5 )M p (4) x 

u a (5)Y(v 4 + a 4 75)M 6 ) = V{XP,pa,ar) A(A, p, cr; 1, 2, 3, 4, 5, 6) 

a,P,T=± 

where 

V(\(3, pa, or) = P X pP pa P aT VfV 2 a V 3 T (8) 

= vi ± a,\ 
V 2 = (v 2 ± a 2 )(v 3 ± a 3 ) 
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= Vi ± 04 
A(A,p,a;l,2,3,4,5,6) = 
m a (1) 7 ^a(2) x 

m p (3)7 m (^i +fo)luU p (4:) x 
^(5)7^(6) . 

Here we have disregarded the particle/antiparticle distinction since it is al- 
ready implied by the assignment of the external momenta. The helicity 
labels X,p,a = ± determine the helicity of both external legs on a given 
fermion line. Using the Weyl-van der Waerden formalism for helicity am- 
plitudes (or, equivalently, the Dirac formalism of 0), the expression A 
can easily be calculated p. It turns out that, for each permutation of the 
fermion momenta, all helicity combinations can be computed using only four 
independent complex functions. 

The numerator in the non-Abelian diagrams can also be written in terms of 
the function A: 

u x (l)'j a ux(2) ^(3)7^(4) ^(5)7^(6) 
x 2 {g^( Pl + p 2 y + g a »(p 5 + n y + g^(p 3 + p 4 ) a } 
= A(\, p, cr; 1, 2, 3, 4, 5, 6) - A(a,p, A; 5, 6, 3, 4, 1, 2) . (9) 

Thus, for massless fermions, every helicity amplitude consists of a sum of 
very systematic, and relatively compact, expressions. 

When four quarks are present in the final state, one has to add the con- 
comitant QCD production channels, and also the production of a quark pair 
and two gluons, since both types of final states will appear as jets. The 
former contribution, which we call interfering QCD background, is easily 
implemented once all electroweak diagrams have been computed as shown 
before. In fact, it is enough to add gluons wherever photons connect quark 
lines (of course the correct QCD coupling and colour structure should be 
taken into account). Finally, the latter process can be efficiently computed 
using the recursion relations of ref. j8|. Since it does not interfere with the 
other diagrams, we have written a separate event generator to get this con- 
tribution ||. For the sake of brevity we do not describe it here. However, we 
point out that, as for the Monte Carlo integration, it has been built following 
exactly the same strategy used in EXCALIBUR. 
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The problem of the integration over the final fermion momenta can be 



solved using a multichannel approach [10, If /($) denotes the matrix 
element squared and d& the 8-dimensional massless phase space integration 
element, one has to compute 

a = J /($) d$ 9{cuts) , (10) 

where 8 (cuts) stands for any kind of experimental cut, that, in a Monte Carlo 
approach, is implemented by simply putting /($) = in the unwanted region 
of the phase space. 

In order to reduce the variance of the integrand, and therefore the Monte 
Carlo error, it is convenient to introduce an analytically integrable function 
g($), called the local density, that exhibits approximately the same peaking 
behaviour of /($) and is unitary, that is, a normalized probability density: 

g@)d$ = l . (11) 

By multiplying and dividing the integrand by <?(<&), the cross section can be 
rewritten as follows 

a = J w(${p)) dp6(cuts) (12) 
where the new integrand 

w{*(fl) = (13) 

is a smoother function of the new set of variables {p{\ defined by 

dp = g(&)d$ 
< pi < 1 (14) 

so that the variance of w(p) is smaller than the variance of /($)• 
When the peaking structure of the matrix element squared is very rich one 
set of new integration variables {p^} can only describe well a limited number 
of peaks. Therefore a multichannel approach is required in which 

N r - - 

#) = j = l\ ft ($) , 5>* = 1 , j gi ($)d$ = l , (15) 

i=i 
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and where every <7i($) describes a particular peaking structure of /(<&). Note 
that the conditions on the a.i and <7i($) ensure unitarity of the algorithm, 
i.e. probability is explicitly conserved at each step of the algorithm, without 
additional normalization factors at any stage. 

The numbers cti are called a-priori weights and, although their numerical 
values are in principle unimportant, they can be used, in practice, to reduce 
the Monte Carlo error [ [PH| . 

In EXCALIBUR we have dealt with the problem of the construction of the 
in a very modular and systematic way. Firstly, we have singled out all 
possible kinematical diagrams occurring in a four-fermion final state (see fig. 
2). They are pictures, inspired by the Feynman diagrams, which represent 
the various peaking structures of the matrix element and indicate which 
variables are most appropriate to a given The explanation of the 

pictures will be given in section 3.1. Secondly, we have written all building 
blocks (that is subroutines) necessary for the calculation. Finally, we have 
put them together to form the #«($)• 

QED corrections are implemented using the structure-function formalism 
[|l|, [|. Each of the incoming fermions is assumed to have its energy degraded 
by the emission of photons parallel to the beam. For the energy distribution 
of the fermion after radiation we take a structure function $ that incorporates 
the leading log O (a) and O (a 2 ) initial state radiation with exponentiation 
of the remaining soft-photon effects. Its expression can be found in ||. Our 
model for the total radiative cross section is then 

i i 

a(s) = J J dxi dx 2 $(^i) &(x 2 ) a (xix 2 s) (16) 
o o 

where cr is the non-radiative cross section and X\, x 2 represent the energy 
content of the incoming fermions after radiative emission. This provides an 
adequate description of the leading QED effects. 

3 Program Structure 

We shall now describe in some detail the salient features and strategies 
adopted in EXCALIBUR. The Program consists of two parts: the evaluation of 
the matrix element and the event generation. Both steps require an initial- 
ization, according to the chosen final state. Roughly speaking it means that 
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the Feynman diagrams and the kinematical channels have to be built. This 
initialization is done in SUBROUTINE SETPRO, the matrix element is evaluated 
in SUBROUTINE DIAGA and SUBROUTINE MATRIX, while nearly all the rest is 
devoted to the event generation and Monte Carlo integration. 

3.1 Subroutine SETPRO 

We already described the algorithm used to construct the Feynman dia- 
grams through a big do loop over all 36 permutations of the six fermion 
momenta. In SUBROUTINE SETPRO the variable KPERM(1 : 6 , 1 : 36) explicitly 
contains all these permutations, and IPHASE(1:36) the corresponding rela- 
tive phase. The constructed Abelian (non Abelian) diagrams are stored in 
J J ( 1 : 16 , NDAB) ( JN ( 1 : 16 , NNAB) ) , where the first index contains information 
about the particles involved in the process, the vector boson propagators and 
the momenta assignment, and the second one enumerates each diagram. In 
the output each constructed diagram is printed out together with its list num- 
ber. For particular studies or checks, we give the possibility to switch off dia- 
grams. This can be achieved by putting the variables KAO(I)=0 (K10(J)=0), 
for the corresponding unwanted Abelian (non Abelian) diagrams (1=1 :NDAB, 
J=1:NNAB). SUBROUTINE SETPRO also contains the the input parameters of 
the program. They are a (ALPHA), a s (ALS, relevant for 4 quarks final states), 
M z (ZM), M w (WM), sin 2 # w (STH2), Y w (WW) and Y z (WZ). The statistical fac- 
tor STATFAC and the colour factor FCOL are evaluated according to the chosen 
final state. Furthermore, the coupling combinations Vf 1 of eq. |] (and those 
occurring in the non Abelian case) are computed. It may happen that, for a 
particular helicity combination, one or more of the are zero. In the latter 
cases there is no point in computing the corresponding function A (see eq. 
0). As a result, less than four independent complex functions are required to 
evaluate the spinorial part of the diagram. In order to have a fast evaluation 
of the matrix element those cases have to be excluded. This is achieved by in- 
troducing two occupation matrices NC(1:36) and N0C(1 :36, 1 :4). For each 
of the 36 permutation, NC is set zero if the corresponding permutation does 
not give any Feynman diagram, and, if it does, NOC indicates which complex 
functions are needed. Through C0MM0N/AREA3/ these matrices are passed to 
SUBROUTINE DIAGA, where only those helicity combination for which NOC and 
NC are different from zero are computed. 

Two more operations are performed in SUBROUTINE SETPRO, namely the 
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choice of the kinematical channels for the Monte Carlo integration and the 
computation of the QCD interfering background. 

We singled out a maximum number of 26 kinematical channels. They are 
given in fig. 2, together with the name of the corresponding subroutines in 
EXCALIBUR, and are inspired by all possible occurring Feynman diagrams. 
Fermionic lines have an arrow, a wavy line represents a photon and a dashed 
line can be either Z or W (this gives 26 channels). Solid lines connect 
topological equivalent points. That is a t-channel solid line means isotropic 
angular distributions between the connected fermions while a s- channel solid 
line stands for photon or massive vector boson propagators. Since they only 
give rise to an s dependent behaviour, the peaking structure relevant for the 
integration over the final momenta is not affected by them. As an example, 
with those conventions it is easy to recognize that the last channel RAMB04 
represents an isotropic 4 body decay. If the inspiring Feynman diagrams ex- 
ist, the variables NCHA ( 1 : 26 , 1 : 48) are set equal to one, where the first index 
runs over the possible channels and the second one labels the permutation 
of the final momenta. The number 48 is explained as follows. There are 
24 permutations of the four final momenta but, for some channel, the case 
where the initial state labels 1 and 2 refer to e + and e~ respectively must be 
distinguished from the case where they refer to e~ and e + . This gives 48 pos- 
sible permutations. Depending on the topology, there are symmetries among 
the final momenta that have to be taken into account in order to have a 
minimum number of kinematical channels. For example, in channel N0NAB1, 
permutation 3456 of the final momenta is equivalent to permutation 5634. 
This symmetrization is automatically performed by the program. When the 
initialization in SUBROUTINE SETPRO is completed, variable NCT indicates the 
number of found kinematical channels. In the output file they are printed 
out together with a list number I. An array has been introduced (NCTO), such 
that putting NCTO (I) = excludes channel with number I (1=1: NCT) when 
the Monte Carlo integration is performed. This can be used to increase the 
speed of the program, by switching off those channels for which the procedure 
of self-optimization (see below) gives very small a-priori weights. 
The interfering QCD background is added as an extra contribution propor- 
tional to the ratio -^qi (the variable GRAP), for those amplitudes where a 
photon connects two quarks of charge Q and Q'. 

This concludes the description of SUBROUTINE SETPRO. Since all possible ini- 
tializations are performed there, the structure of the rest of the program can 
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be simple and fast. 



3.2 Subroutines DIAGA and MATRIX 

In SUBROUTINE DIAGA numerators and denominators of all found Feynman 
diagrams for which NOC and NC are non vanishing are computed at once using 
the Weyl-van der Waerden formalism. It means that, for each generated 
event, SUBROUTINE DIAGA is called just once and not n times, where n is the 
total number of Feynman diagrams. As for the computational speed, this is 
very important. 

In SUBROUTINE MATRIX (SQUAREM) the matrix element squared (SQUAREM) is 
calculated by putting together the numerators and denominators computed 
in SUBROUTINE DIAGA and the coupling combinations of eq. |8] evaluated in 
SUBROUTINE SETPRO. Since computing the colour factor and QCD interfering 
background in a four quark final state with colour labels i, j, I and m requires 
the part of the amplitude proportional to 5ij5i m to be distinguished from that 
proportional to SaSj m [|J, the constructed amplitudes in SUBROUTINE MATRIX 
take care of both contributions separately. 



3.3 Phase space generation and integration 

In the MAIN of EXCALIBUR the variables XR1 and XR2, representing the energy 
content X\ and X2 of the incoming fermions after radiative emission (eq. [T^), 
are generated. Then, the initial configuration of the momenta in the center 
of mass frame of the event after ISR is set calling SUBROUTINE MOMSET and 
the cuts imposed on the momenta in the Lab frame are rewritten in terms 
of cuts in the center of mass frame. The kinematical channels are called 
using SUBROUTINE ADDRESS (LFLAG , NC , NN , D J) . When LFLAG is set 0, the 
channel number NC, with the momenta permutation labelled by NN, is used 
for generating the momenta and computing the local density DJ. If LFLAG= 
1 the actual momenta configuration is used to compute DJ. The choice of the 
channel to use is performed, as in ref . JTO] , on the basis of the actual values of 



the a-priori weights by defining the cumulative numbers = ax + ... + 0$, 
taking a random number uniformly distributed between and 1 and choosing 
channel i if < z < /3j. 

In SUBROUTINE MOMARRAY the generated four-momenta are put in a big array 
PM(0:4, 0:900) and stored in 
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COMMON/MOMENTA/ROOTS , XR1 , XR2 , PM (0 : 4 , : 900) 



(ROOTS is the center of mass energy of the event). The first index refers to 
the component of the momenta (0 represents the energy and 4 is the four 
momentum squared). As for the second index, the following self-explanatory 
conventions are used: 

PM(I,34) = PM(I,43)= PM(I,3)+PM(I,4) etc. 

PM(I,643) = PM(I,346) = ... = PM(I ,6)+PM(I ,4)+PM(I ,3) etc. 

Besides, but only if the first digit refers to an incoming momentum (notice 
the correspondence 7 — > 1, 8 — > 2) 

PM(I,734) = PM(I,1)-PM(I,3)-PM(I,4) etc. 
PM(I,851) = PM(I,2)-PM(I,5)-PM(I,1) etc. 

Each channel is constructed in a very modular way by putting together basic 
subroutines that describe different parts of its peaking structure. In ref. M an 
example of the construction of channels BREMB2 and C0NVER2 of fig. 2 is given. 
There are 10 of these basic subroutines. They are the building blocks of the 
whole generation procedure. For the sake of brevity we do not list them here. 
They are well commented in the program. We only notice that, in building 
the kinematical channels, every t-channel exchanged massive vector boson is 
always assumed to give a flat angular distribution between the initial and the 
final fermion. This is done in order to avoid proliferation in the number of 
channels. In our experience, this gives a very good approximation at center 
of mass energies up to 500 GeV, a good approximation at higher energies up 
to 1 TeV and may cause large Monte Carlo errors at 2 Tev. Of course the 
Monte Carlo program remains correct, but higher statistics runs are required. 
However, adding channels to map this high energy kinematical behaviour is 
trivial, because EXCALIBUR already contains all needed ingredients. 

As far as the self-optimization of the integration is concerned, a detailed 
description of the iterative algorithm implemented in EXCALIBUR may be 
found in ref. [JTTJ] . Here we point out that two variables have to be chosen by 
the user, namely the maximum number of iterations ISTEPMAX (in the input 
list) and the number of point N0PT used for the self-optimization (in the MAIN 
of the program). Then, for each iteration, N0PT/ISTEPMAX points (including 
zero-weight events) are used to compute the a-priori weights. We found 
that, with 4-5 hundred thousand points, a good choice is N0PT= 100,000 and 
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ISTEPMAX= 10. However, when very stringent cuts are applied, the majority 
of the events falls outside the allowed region, so that the ratio NOPT/ISTEPMAX 
may be a very small number. This causes a bad estimate of the best a-priori 
weights to be used. In those cases it is convenient to either increase NOPT or 
decrease ISTEPMAX. 

In the input list one has to specify the set of standard cuts as specified 
in the next section. Any other type of cut must be implemented directly in 
SUBROUTINE CUTS (LNOT) , where 

C0MM0N/AREA10/PM1 (0 : 4 , 1 : 6) , PM4 ( 12 : 65) , 0MCT1 (1:6,3:6) 

contains the four momenta computed in the Lab frame (PM1), the invariant 
mass squared among all possible particles pairs (PM4) and the quantities 
1 — cos% (0MCT1). If the event is rejected LN0T= 1, and the weight is put 
equal to zero. 

Finally, all weights (computed as in eq. |13|) are collected using SUBROUTINE 
INB00K and the Monte Carlo results called through SUBROUTINE 0UTB0K. 

4 Input 

The meaning of the input parameters is the following: 
NPROCESS (INTEGER) 

The number of processes to be computed. 
N (INTEGER) 

The number of points for the Monte Carlo integration. 
ISTEPMAX (INTEGER) 

The number of iterations for optimizing the a-priori weights. 

OUTPUTNAME ( CHARACTER* 1 5 ) 

The name of the output file. 

KREL( INTEGER) 
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It selects the signals. If KREL= all Feynman diagrams are taken into account. 
If KREL= 1-5 only those leading to reactions of eqs. (HD - ©- 

LQED( INTEGER) 

It includes (1) or excludes (0) ISR. 
R00TSMUL (REAL*8) 

The total energy of the colliding e + and e~. All energies are in GeV. 
SHCUT(REAL*8) 

Minimum value of the invariant mass squared of the event after QED radia- 
tion. 

ECUT(3) (REAL*8) 

Minimum energy of particle number 3. 
ECUT(4) (REAL*8) 

Minimum energy of particle number 4. 
ECUT(5) (REAL*8) 

Minimum energy of particle number 5. 
ECUT(6) (REAL*8) 

Minimum energy of particle number 6. 
SCUT (3, 4) (REAL*8) 

Minimum value of (j>(3) + p(4)) 2 . All invariant masses are in GeV 2 . 

SCUT (3, 5) (REAL*8) 

Minimum value of (p(3) +p(5)) 2 . 

SCUT (3, 6) (REAL*8) 

Minimum value of (p(3) + p(6)) 2 . 

SCUT (4, 5) (REAL*8) 
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Minimum value of (p(4) +p(5)) 2 . 

SCUT (4, 6) (REAL*8) 

Minimum value of (p(4) +p(6)) 2 . 

SCUT (5, 6) (REAL*8) 

Minimum value of (p(5) +p(6)) 2 . 

CMAX(1,3) (REAL*8) 

Maximum value of cos 6* between particle 1 and 3. 
CMAX(1,4) (REAL*8) 

Maximum value of cos 9 between particle 1 and 4. 
CMAX(1,5) (REAL*8) 

Maximum value of cos 9 between particle 1 and 5. 
CMAXC1.6) (REAL*8) 

Maximum value of cos 9 between particle 1 and 6. 
CMAX(2,3) (REAL*8) 

Maximum value of cos 9 between particle 2 and 3. 
CMAX(2,4) (REAL*8) 

Maximum value of cos 9 between particle 2 and 4. 
CMAX(2,5) (REAL*8) 

Maximum value of cos 9 between particle 2 and 5. 
CMAX(2,6) (REAL*8) 

Maximum value of cos 9 between particle 2 and 6. 
CMAX(3,4) (REAL*8) 

Maximum value of cos 9 between particle 3 and 4. 
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CMAX(3,5) (REAL*8) 

Maximum value of cos 9 between particle 3 and 5. 
CMAX(3,6) (REAL*8) 

Maximum value of cos 9 between particle 3 and 6. 
CMAX(4,5) (REAL*8) 

Maximum value of cos 9 between particle 4 and 5. 
CMAX(4,6) (REAL*8) 

Maximum value of cos 9 between particle 4 and 6. 
CMAX(5,6) (REAL*8) 

Maximum value of cos 9 between particle 5 and 6. 
PAR (3) (CHARCTER*8) 

Produced fermion with label 3 (to be chosen among ' EL ' , ' NE ' , ' MU ' , ' NM ' , 
'TA', 'NT', 'DQ', 'UQ', 'SQ', 'CQ', 'BQ', ' TQ ' ) . 

PAR (4) (CHARCTER*8) 

Produced antifermion with label 4. 

PAR (5) (CHARCTER*8) 

Produced fermion with label 5. 

PAR (6) (CHARCTER*8) 

Produced antifermion with label 6. 

5 Test Run Output 

To conclude our description, we give an example of a typical calculation 
that can be performed with EXCALIBUR. One should be able to reproduce 
this output within the estimated Monte Carlo error (small differences may 
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occur because the quasi-random number generator used in the program is 
not strictly portable). Using an input file as follows 



1 number of energy points 

250000 number of Monte Carlo points 

10 number of iterations in a. p. weights optimization 

output output program name 

krel (signal: 0,1,2,3,4,5) 

1 lqed (0 or 1) 
190. dO total energy (GeV) 

O.dO cut on reduced inv. mass squared after ISR 

O.dO ecut_3 

O.dO ecut_4 

20. dO ecut_5 

20. dO ecut_6 

O.dO scut_34 

O.dO scut_35 

O.dO scut_36 

O.dO scut_45 

0. dO scut_46 
100. dO scut_56 

1. dO cmax_13 
l.dO cmax_14 
0.9d0 cmax_15 

0. 9d0 cmax_16 

1 . dO cmax_23 
1 . dO cmax_24 
. 9d0 cmax_25 

. 9d0 cmax_26 

1 . dO cmax_34 
1 . dO cmax_35 
1 . dO cmax_36 
1 . dO cmax_45 
1 . dO cmax_46 
. 9d0 cmax_56 

mu produced fermion (3) 

nm produced antifermion (4) 
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uq 
dq 



produced fermion (5) 
produced antifermion (6) 



and the values ALPHA= l./128.,ZM= 91.16, WM= 80.22, STH2= 0.226 , WW= 
2 . 03, WZ= 2 . 53 we get the following output file 

output 

All Feynman diagrams 

sqrt(s) = .190000D+03 
n.points = 250000 
istepmax = 10 

energy cuts with ecut_3 = .0 

ecut_4 = .0 

ecut_5 =20.0 

ecut_6 =20.0 

cut on s*xlr*x2r = .0 

mass cuts with scut_34 = .0 

scut_35 = .0 

scut_36 = .0 

scut_45 = .0 

scut_46 = .0 

scut_56 = 100.0 



angle cuts with 



cmax_ 


.13 




1.0 


cmax_ 


.14 




1.0 


cmax_ 


.15 




.9 


cmax_ 


.16 




.9 


cmax_ 


.23 




1.0 


cmax_ 


.24 




1.0 


cmax_ 


.25 




.9 


cmax_ 


.26 




.9 


cmax_ 


.34 




1.0 
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I.S.R. INCLUDED 



cmax_35 = 1.0 

cmax_36 = 1.0 

cmax_45 = 1.0 

cmax_46 =1.0 

cmax_56 = . 9 



s~2_thet 

Z-mass 

Z-width 

W-mass 

W-width 

1/ alpha 

alpha_s 



.226000D+00 
.911600D+02 
.253000D+01 
.802200D+02 
.203000D+01 
. 128000D+03 
. 103000D+00 



process: antiel(l) el(2) > mu(3) antinm(4) uq(5) antidq(6) 



abelian diagrams 



[el(l) 
[el(l) 
[el(l) 
[el(l) 
[mu(3) 
[mu(3) 
Euq(5) 
Euq(5) 



,el(2)] 
,el(2)] 
,el(2)] 
,el(2)] 
,nm(4)] 
,nm(4)] 
,dq(6)] 
,dq(6)] 



[mu(3) 
[mu(3) 
[uq(5) 
[uq(5) 
[uq(5) 
[uq(5) 
[el(l) 
[mu(3) 



,mu,nm(4)] 
,mu,nm(4)] 
,uq,dq(6)] 
,uq,dq(6)] 
,dq,dq(6)] 
,dq,dq(6)] 
,ne,el(2)] 
,nm,nm(4)] 



W [uq(5) 

W [uq(5) 

W [mu(3) 

W [mu(3) 

Z [el(l) 

G [el(l) 

W [mu(3) 

Z [el(l) 



,dq(6)] 
,dq(6)] 
,nm(4)] 
,nm(4)] 
,el(2)] 
,el(2)] 
,nm(4)] 
,el(2)] 



phase 

ph= 
ph= 
ph= 
ph= 
ph= 
ph= 
ph= 
ph= 



non abelian diagrams 



phase 



1: [uq(5) ,dq(6)] [el (1) , el (2)] [mu(3) ,nm(4)] (WZW) ph= 1 
2: [uq(5) ,dq(6)] [el (1) , el (2)] [mu(3) ,nm(4)] (WGW) ph= 1 

kinematical diagrams 



channel 



permutation 
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annihi2(wm) 

annihi2(wm) 

annihi2(wm) 

annihi2(wm) 

conver3(wm) 

nonabl (wm) 

rambo4 



2 
2 
2 
2 
2 
2 
2 



********** weights analysis ********** 



*** variable number 1 ************ 



sum(w**0) 

sum(w**2) 

sum(w**4) 

maximum 

no . weights=0 

estimator x: 

estimator y: 

estimator z: 

average estimate 



.250000D+06, 
. 254349D+06 , 
. 43441 1D+07 
.221013D+02, 
59062, 

. 543028D+00 

.289008D-05 

.730267D-15 

: .543028D+00 



sum(w**l) 
sum(w**3) 

max. in buffer 
no . weights<0 



+\- .170002D-02 
variance estimate: .289008D-05 

+\- .270235D-07 
efficiency for all weights : 2.457 % 

efficiency for nonzero weights : 3.217 % 
overshoot factor of histogram : 1.487 
the distribution of the nonzero weights: 
50, log scale; entries under, inside, over: 



. 135757D+06 
.803758D+06 

. 148670D+02 




190928 



10 



1487E+01 


.1631E+06 




2973E+01 


.2157E+05 




i 


4460E+01 


.4936E+04 




i 


5947E+01 


. 1089E+04 




i 


7433E+01 


. 1470E+03 




i 


8920E+01 


.3600E+02 




i 


1041E+02 


. 1600E+02 


i******* 


i 
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. 1189E+02 .1100E+02 i****** 
. 1338E+02 .7000E+01 i***** 
. 1487E+02 .5000E+01 i***** 



i 
i 
i 



differences in the computation 
of the a-priori weights: 



diff ( 


1) = 


1.95407687917905 


diff ( 


2) = 


.891414805616945 


diff ( 


3) = 


.7677390572607608 


diff ( 


4) = 


.6967978248575471 


diff ( 


5) = 


.5755243287941117 


diff ( 


6) = 


.7280986181406474 


diff ( 


7) = 


.5755417220460467 


diff ( 


8) = 


.5925286045906328 


diff ( 


9) = 


.5632310939026423 


diff ( 


10) = 


.598279075660049 


diff ( 


11) = 


.5921705560137981 


diff ( 


12) = 


.5181929881140225 



a-priori weights: 



1 : .899493D-03 

2 : .165631D-03 

3 : .114732D-03 

4 : .460900D-03 

5 : .881137D+00 

6 : .117221D+00 

7 : .817789D-06 



After information about input parameters and imposed cuts, the program 
prints out the used Feynman diagrams and kinematical channels. Then, the 
analysis of the weights giving the Monte Carlo estimate of the cross section 
(variable number 1) follows. In particular various sums of the weights to 
powers 0-4 are given as well as the maximum weight and that one in the 
buffer (that is in the interval of values used in the histogram that shows the 
weight distribution). The quantity x is the estimator of the average of the 
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distribution defined, for N weights Wi, as 



N 



y is the estimator of the variance 



1 

N 



N 



N(N - 1) 



(17) 



(18) 



and z is an estimator for the variance of the variance, so that the error on the 
average and variance estimates are ^fy and yfz respectively. In the example 

and the 



<w> 



a = .5430 ± .00170 pb. As usual the efficiency is defined as 

L J max(w) 

overshoot factor is the ratio between the maximum weight and the maximum 
weight written in the buffer. In the histogram the 190928 non zero weights 
are displayed according to their abundance in bins. Finally, the variables 
D of ref. [J l|| (that measure, at each step in the optimization procedure, 
how well the actual set of a-priori weights approximates the behaviour of 
the optimal set) are printed out, together with the found best set of a-priori 
weights. 
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Figure 2: kinematical diagrams in EXCALIBUR. 
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